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Minimum time control of the rocket attitude reorientation 
associated with orbit dynamics 

Jiamin Zhu* Emmanuel Trelat^ Max Cerfl 


Abstract 

In this paper, we investigate the minimal time problem for the guidance of a rocket, 
whose motion is described by its attitude kinematics and dynamics but also by its orbit 
dynamics. Our approach is based on a refined geometric study of the extremals coming from 
the application of the Pontryagin maximum principle. Our analysis reveals the existence of 
singular arcs of higher-order in the optimal synthesis, causing the occurrence of a chattering 
phenomenon, i.e., of an infinite number of switchings when trying to connect bang arcs with 
a singular arc. 

We establish a general result for bi-input control-affine systems, providing sufficient condi¬ 
tions under which the chattering phenomenon occurs. We show how this result can be applied 
to the problem of the guidance of the rocket. Based on this preliminary theoretical analy¬ 
sis, we implement efficient direct and indirect numerical methods, combined with numerical 
continuation, in order to compute numerically the optimal solutions of the problem. 

Keywords: Coupled attitude orbit problem; optimal control; Pontryagin maximum principle; 
shooting method; continuation; chattering arcs. 
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1 Introduction 

The optimal control of orbit transfer (see, e.g., [7, 9, 26]) and attitude reorientation (see, e.g., 
[4, 31, 32]) for spacecrafts have been extensively studied in the past few decades. The optimal 
control problem of orhit transfer focuses mostly on how to move the spacecraft from one orbit 
or point to another orbit or point by using minimum energy, while the optimal control problem 
of attitude reorientation is mainly devoted to determine how to change the pointing direction of 
the spacecraft in minimum time. In the existing literature, these two optimal control problems 
are considered separately in general. From the engineering point of view, for most satellites, it 
is appropriate to design separately the control laws for the orbit movement and for the attitude 
movement. However, for the rockets, the trajectory is controlled by its attitude angles: the way 
to make the rocket follow its nominal trajectory is to change its attitude angles, and therefore it 
is desirable to be able to determine the optimal control subject to the coupled dynamical system. 
Though the control of the coupled problem was also studied in many previous works (see, e.g., 
[23, 17, 20]), it does not seem that the problem has been investigated in the optimal control 
framework so far. 

In this paper, we consider the time minimum control of the attitude reorientation coupled with 
the orbit dynamics of a rocket, denoted in short (MTCP). The chattering phenomenon that may 
occur according to the terminal conditions under consideration, makes in particular the problem 
quite difficult. Chattering means that the control switches an infinite number of times over a 
compact time interval. Such a phenomenon typically occurs when trying to connect bang arcs 
with a higher-order singular arc (see, e.g., [13, 24, 35, 36]). In [36], we studied the planar version 
of (MTCP), where the system consists of a single-input control-affine system, and we established 
as well the occurence of a chattering phenomenon and that the chattering extremals are locally 
optimal in (7° topology.* 

A second important difficulty in (MTCP) is due to the coupling of the attitude movement 
with the orbit dynamics. Indeed the system contains both slow (orbit) and fast (attitude) dy- 

*A trajectory x(-) is said to be locally optimal in topology if, for every neighborhood V of x(-) in the state 
space, for every real number rj so that |t 7 | e, for every trajectory x(-), associated to a control v on [0, T + r;], 
contained in W, and satisfying xfO) = SfO) = xq, x{T + rj) = :r(T), there holds C{T + rj^r) ^ C(T, u), where C is 
the cost functional to be minimized. 


2 











namics. This observation will be particularly important in order to design appropriate numerical 
approaches. 

In order to analyze the extremals of the problem, we use geometric optimal control theory (see 
[1, 30, 33]). The Pontryagin maximum principle and the geometric optimal control, especially 
the concept of Lie bracket, will be used in this paper in order to establish an existence result of 
the chattering phenomenon. More precisely, based on the Goh and generalized Legendre-Clebsch 
conditions, we prove that there exist optimal chattering arcs when trying to connect a regular arc 
with a singular arc of order two. 

There exist various numerical approaches to solve an optimal control problem. The direct 
methods (see, e.g., [3]) consist of discretizing the state and the control and thus of reducing the 
problem to a nonlinear optimization problem (nonlinear programming) with constraints. The in¬ 
direct methods consist of numerically solving a boundary value problem obtained by applying 
Pontryagin maximum principle (PMP, see [29]), by means of a shooting method. There exist also 
mixed methods that discretize the PMP necessary conditions and use then a large-scale optimiza¬ 
tion solver (see, e.g., [2]). Since these numerical approaches are not easy to initialize successfully, 
it is required them to combine with other theoretical or numerical approaches (see the survey [33]). 
Here, we will use numerical continuation, which has proved to be very powerful tool to be combined 
with the PMP. For example, in [10, 15, 25], the continuation method is used to solve difficult orbit 
transfer problems. 

However, due to the chattering phenomenon, numerical continuation combined with shooting 
cannot give an optimal solution to the problem for certain terminal conditions for which the optimal 
trajectory contains a singular arc of higher-order. In that case, we propose sub-optimal strategies 
by using direct methods computing approximate piecewise constant controls. It is noticeable that 
our indirect approach can also be adapted to generate sub-optimal solutions, by stopping the 
continuation procedure before its failure due to chattering. This approach happens to be faster 
than the direct approach, and appears as an interesting alternative for practice. 

From the engineer point of view, the theoretical analysis as well as the numerical strategies and 
the way to design them (in particular, the design of the problem of order zero) are strongly based 
on the fact that the orbit movement is much slower than the attitude movement. 

The paper is organized as follows. In Section 2, we describe the mathematical model of the 
system consisting of the attitude dynamics, of the attitude kinematics, and of the orbit dynamics. 
In Section 3, we recall the Pontryagin maximum principle and some higher necessary conditions of 
optimality (Goh and generalized Legendre-Glebsch conditions) for bi-input control affine systems. 
Based on these necessary conditions of optimality, we establish a result on the existence of the 
optimal chattering extremals. In Section 4, we analyze the regular and singular extremals of the 
problem (MTCP). For the regular extremals, we classify the switching points and state some 
useful properties. For the singular extremals (which are of order two), we show that the chattering 
phenomenon occurs for the problem (MTCP) by using the results given in the previous section. 
In Section 5, we propose a numerical approach to solve the problem (MTCP) by implementing 
numerical continuation combined with shooting. Numerical results are given in Section 6. 

2 Model and problem statement 

The problem is to control the attitude movement coupled with the orbit dynamics in the 
launching ascent stage for a rocket. In this paper, we take the system parameters of the rocket 
Ariane 5. In order to keep the stability of the rocket along the flight, the attitude maneuver should 
be moderate, i.e., at most ±20 degrees, and then it is possible to use Euler angles to model the 
attitude of the engine. In this section, we first define the coordinates systems, and then we give 
the equations of the attitude dynamics, of the attitude kinematics and of the orbit dynamics. The 
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model consists of eight ordinary differential equations: three for the components of the velocity 
vector, three for the Euler angles and two for the components of the angular velocity vector. 


2.1 Coordinate systems 


Throughout the paper, we make the following assumptions: 

• The Earth is a sphere and is fixed in the inertial space, i.e., the angular velocity of the Earth 
is zero, which means that Wei = 0- 

• The position of the rocket remains the same during the maneuver of the rocket. 

• The rocket is an axial symmetric cylinder. 

• The aero-dynamical forces are zero. 

• The rocket engine cannot be shut off during the flight and the module of the thrust force is 
constant, taking its maximum value, i.e., T = T^ax- 

The unit single-axis rotation maps i?i(cr): M —>■ for ct £ R, i = x, y, z are defined by 


/I 

0 

0 \ 

/ COSfJ 

0 

— sin a\ 

Rx{<j) = 0 

cos a 

sin a 

1 ) Rvi<^) = 0 

1 

0 

\o 

— sina 

cosa / 

\sincr 

0 

COSCr / 


( cos a sin a 
— sin a cos a 
0 0 



Eor a given vector e £ taking Ri{a)e means to rotate the vector e with respect to the axis i 
by an angle of a. With this definition, we next introduce the coordinate frames that will be used 
throughout the paper. 

The Earth frame Sg = {Xg,yg,Zg) is fixed around the center of the Earth O. The axis Zg 
points to the North pole, and the axis Xg is in the equatorial plan of the Earth pointing to the 
equinox. 

The launch frame Sr = {xR,yR, zr) is fixed around the launch point Or (where the rocket 
is launched). The axis xr is normal to the local tangent plane, pointing to the launch direction 
(here we assume that the rocket is vertically launched, i.e., the launch direction is perpendicular 
with the local tangent plane), and the axis zr points to the North. As shown in Eigure 1 (a), the 
launch frame is derived from the Earth frame by two ordered unit single-axis rotations Rz{(-r) and 
Ryi-LR), 

„ Rz(f-R) Ryi-La) 

Sg -)> O -)> Sr 

where iR and Lr are the longitude and latitude of the launch point, respectively. 

The body frame Sb = {xb,yb,Zb) is defined as follows. The origin of the frame Ob is fixed 
around the mass center of the rocket, the axis Zb is along the axis-symmetric axis of the rocket, and 
the axis Xb is in the cross-section. The body frame can be derived by three ordered unit single-axis 
rotations from the launch frame, as shown in Figure 1 (b), 


Sr 


Ryi^) Rxi^l^) Rzi4*) Q 
- > o - > o - y Sb 


where 9 is the pitch angle, "0 is the yaw angle and 0 is the roll angle. Therefore, the transformation 
matrix from Sr to Sb is 


LbR =Rz{4>)Rx{tp)Ry{d) 


( cos 9 cos 0 -I- sin 9 sin 0 sin 0 
— cos 9 sin 0 -|- sin 9 sin 0 cos 0 
sin 9 cos 0 


cos 0 sin 0 
cos 0 cos 0 
— sin0 


— sin 9 cos 0 -I- cos 9 sin 0 sin 0^ 
sin 9 sin 0 -|- cos 9 sin 0 cos 0 
cos 9 cos 0 


( 1 ) 
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Figure 1: Coordinate systems and relations. 


and the transformation matrix from Sf, to is Lm, = 

2.2 Attitude dynamic equations 

The attitude dynamics are written in vectorial form in the body frame Sb as 

{Iuj)b = —{^)b A {I(2)b + iM)b, (2) 

where / is the inertia matrix, lo is the absolute angular velocity vector, i.e., the angular velocity 
of the rocket with respect to the inertial space, and M is the control torques introduced by the 
rocket thrust. The index {■)b means that the vectors are expressed in the body frame Sb- 
Setting {I)b = diag(4,4,/^), (w);, = (w^;,Wy,and {M)b = , (2) gives 

{ ^x^x — Iz)^y^z T 

~ (A Ix']^x^z T Tfy, (3) 

Iz^z = {Ix — Iy)tOx^y + 

The control torque M is the cross product of the thrust vector T and of its moment arm L. The 
moment arm is the vector from the center of mass Ob to the force acting point Op, given here by 
{L)b = (0, 0, —l)^- Moreover, as shown in Figure 2 (a), the thrust force vector is 

T = (—T sin/i cos C, —T sin/i sin C, T cos/i)^, 

where T = T^ax, M G [0,/imax], and ( G [—7r,7r]. The control torque is then 

{M)b = {L)b A {T)b = {—Tlsm^sm(, Tls'mficos(, 0)^. 

By assumption, the rocket is axial symmetric, and hence A = Iv- Assume that WzlO) = 0, and let 
b = T-maxl/Ix- Then (3) gives 

Wx = —fesin^sin^, uiy = b sin fi cos C., 

with uiz = 0. 

According to the parameters of the rocket engine, Umax is less than 10 degrees and thus the 
error between sin^ and ^ is less than 0.5%. Therefore, in the model we make the approximation 
sinp, ~ ^ and we define ui = p-cosC, and U 2 = p-sinC, and p = ^max with h = b^max- Hence 

uJx = -bu 2 , ojy = bui- ( 4 ) 
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Figure 2: Thrust in the body Frame 


2.3 Attitude kinematics equations 

Since w is the angular velocity vector of the rocket with respect to the inertial space, it is 
equal to the sum of the angular velocity ujbg of the rocket with respect to the Earth frame, of the 
angular velocity ujge of the Earth frame with respect to the Earth, and of the angular velocity ujei 
of the Earth with respect to the inertial space. According to the assumptions and definitions of 
the frames, it is easy to see that the last two terms are zero, and thus uj = uJbg- Therefore, based 
on the definition of the body frame, the relationship between angular velocity and Euler angles are 

/ 0 \ /tjj cos (l)\ / 0 \ 

= LbR f 0 1 + f -•!/;sin(^ 1 + f 0 1 , 

where L^r is given by (1). Then the equations of the attitude kinematics are 

9 = (ujx sin (j) + uiy cos (j))/cosip, ip = uixcoscp — WySincp, <p = inn ip {ujx sin p + lUy cos (p). (5) 

Therefore, the two equations of (4) and the three equations of (5) describe the attitude movement. 

Note that when ip = 7r/2 + fc7r, fc G N, the Euler angles defined above are not well defined (usual 
singularities of the Euler angles). We assume in this paper that the maneuvers are small enough, 
so that these singularities will not be encountered. 

2.4 Orbit dynamics equations 

The equation of the orbit dynamics in vectorial form is 

= i9)R + + (cd)fl A (E)k - 2(^,, A V)r - (die. A {UJR A r))«, (6) 

at m 

where the notation [pR means that the vector is expressed in the launch frame Sr. The vector 
V is the velocity of the rocket with respect to the Earth frame Sg, and its components in the 
launch frame are Vx, Vy and Vz- The vector {g)R = [gx, Qy, Qz)^ can be approximated by {g)R k, 
(— 50 iO, 0)^, where go is a real number representing the standard gravity (go = 9.8). 
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Due to the fact that the control angle ^ is very small in practice (physical constraints imposed 
by the rocket engine), we assume that the thrust force is along the body axial symmetric axis. 
According to the previous assumptions, the equation of the orbit dynamics (6) becomes 

Ua; = a sin 0 cos'0 + , Vy =—asintp + gy, = acosO cosip + g^, (7) 

where a = Tmax/’^n is constant. Note that this additional assumption is made also because the 
attitude of the rocket is controlled by only a part of the rocket engines, and so the total thrust 
remains almost parallel to the rocket symmetric axis. 

2.5 Minimum time control problem (MTCP) 

Model. The system (4)-(5)-(7) has two control inputs ui and U 2 , and we obtain the system 
Vx = asm0 cosip + gx, Vy =—asmip + gy, = acos6 cosip + gz, 

9 = (ujx sin cp + LVy cos (p) / cosip, ip = ujx cos (p — ujy sin cp, p = (ojx sin p + ujy cos p) tan ip, 

OJx = —bU2, djy = bui. 

( 8 ) 

Defining the state variable x = {vx,Vy,Vz,0,ip,p,uix,u>y), we write the system (8) as the bi-input 
control-affine system 

X = f{x)+uigi{x)+U 2 g 2 {.x), (9) 

where the controls ui and U 2 satisfy the constraint u\ + U 2 ^ 1, and the vector fields /, gi and g 2 
are defined by 


d d d 

f = (asin0 cos0 -I- gx)-K -f {—asinip + gy)-^ -h {a cos 9 cosip -|- gz )-^— 

OVx 


dVy 


dvz 


0 0 0 

-I- {ojx sin0 -I- ojy cosp)/cosip— + {uix cosp — ujy sin0)— -|- tax\p{ujx sin0 -|- ujy cosp) — , 

00 O^jj 0(j) 

gi = b^, g2 = -b^. (10) 


dUJy 


duJx 


Terminal conditions and system parameters. Let Vyoi Vzo, 'Po, Po, xIxq, ujy„, 9f, ipf, 
pf, ujxf and U!yj be real numbers. The initial conditions are fixed to 

l^a:(0) = Vxo, 'yy(O) = Vyg, ?;,^(0) = Vza, 

0 ( 0 ) = 00 , 0 ( 0 ) = 00 , 0 ( 0 ) = 00 , uJx{0)=iXxo, ujy{0)=u;yg. 

The desired final velocity is required to be parallel to the body axis Zb, according to {V{tf))R A 
[zb{tf))ii = 0, and therefore, the constraints on the final conditions are 


Vzf sinipf + Vyj cos9f cosipf = 0, Vzf sm9f — Vxf cos9f = 0, 

9{tf)=9f, lP{tf)=1pf, P{tf)=Pf, UJx{tf)=U:xj, UJy{tf) =UJyj. 

Note that the parallel condition on the final velocity is due to the fact that most rockets are planned 
to maintain a zero angle of attack along the flight. The angle of flight, when the air wind is set to 
zero, is defined as the angle between the velocity and the rocket body axis. 
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Minimum Time Control Problem (MTCP). We set a:o = [vxo.Vy^^Vza.OQ, tpo, (^o, Wyp) G 
M®, and we define the target set (submanifold of K®) 

Ml ={{vx,Vy,Vz,0,ilj,(j),uJx,u}y) I i;^sin'!/'/+ Uycos0/cos'0/= 0, 

Uz sin^/ + ?;j, COS0/ cos'0/= 0, d = df, = 4' = 4'f^ = ^y=^vf}- 

The minimum time control problem (MTCP) consists of steering the bi-input control-affine system 
(9) from x(0) = Xq to the final target Mi in minimum time tf, with controls satisfying the constraint 
itf -I- ^ 1. 

3 Some general results for bi-input control-affine systems 

In this section, we focus on the chattering phenomenon for bi-input control-affine systems with 
control constraints and with commuting controlled vector fields. The results that we are going to 
give are general and will be used in the next section to analyze the problem (MTCP). 

We consider the following general framework. Let M be a smooth manifold of dimension n, let 
xo G M he arbitrary, and let Mi be a submanifold of M. We consider on M the minimal time 
control problem 


{ min t f , 

i(t) = f{x(t))+ui{t)gi{x{t))+U 2 {t)g 2 {x{t)), u={ui,U 2 ) 

\\u(t)\\‘^ = Ui(t)‘^ + U2{tf < 1 , 

x(0) = xo, x{tf) G Ml, tf ^ 0 free, 

where /, gi and g 2 are smooth vector fields on M. 

According to classical results (see, e.g., [11, 33]), there exists at least one optimal solution 
(x(-),u(-)), defined on [0,t/]. 


3.1 Application of the Pontryagin maximum principle 

According to the Pontryagin maximum principle (in short, PMP, see [29]), there must exist an 
absolutely continuous mapping p(-) defined on [0,t/] (called adjoint vector), such that p(t) G 
(cotangent space) for every t G [0,t/], and a real number < 0, with (p(-),p°) ^ 0, such that 

f)ZJ 

x{t) = —{x{t),p{t),p°,u{t)), p{t) = - — {x{t),p{t),p°,u{t)), 


almost everywhere on [0,t/], where H{x,p,p'^,u) = ho(x,p) + uihi{x,p) -I- U 2 h 2 {x,p) + p^ is the 
Hamiltonian of the optimal control problem (13). Here, we have set ho{x,p) = {p, f{x)), hi{x,p) = 
{p,gi{x)), and h 2 {x,p) = {p,g 2 {x)). The maximization condition of the PMP yields, almost 
everywhere on [0,t/], 


L{t) = 


{hi{t),h2{t)) 


m 

ii^wir 


(14) 


whenever ^{t) = {hi{t), h 2 {t)) ^ (0,0). We call $ (as well as its components) the switching 
function. Note that $ is continuous. Here and throughout the paper, we denote by hi{t) = 
hi{x{t),p{t)), with a slight abuse of notation. 

Moreover, we have the transversality condition p(t^) A where Tx(^tf)Mi is the tangent 

space to Ml at the point x(tf), and, the final time tf being free and the system being autonomous, 
we have also ho{x{t),p(t)) + ||$(t)|| +p° = 0, Vt G [0,t/]. 





The quadruple (x(-),p(-),p^,u(-)) is called an extremal lift of x(-). An extremal is said to be 
normal (resp., abnormal) if < 0 (resp., p^ = 0). 

We say that an arc (restriction of an extremal to a subinterval I) is regular if ||$(t)|| 7 ^ 0 along 
I. Otherwise, the arc is said to be singular. Note that a singular extremal may be both normal or 
abnormal. We will see in Section 4.2 that the singular extremals of the problem (MTCP) must 
be normal. 

A switching time is a time t at which <f>(t) = (0,0), that is, both hi and ^2 vanish at time t. 
An arc that is a concatenation of an infinite number of regular arcs is said to be chattering. The 
chattering arc is associated with a chattering control that switches an infinite number of times, 
over a compact time interval. A junction between a regular arc and a singular arc is said to be a 
singular junction. 


3.2 Computation of singular arcs, and necessary conditions for optimal¬ 
ity 

We next define the order of a singular control, since it is important to understand and explain 
the occurence of chattering. This concept is related to the way singular controls are computed, 
and since it is a bit technical to define, we start with a preliminary quite informal discussion. Here 
and throughout the paper, we use the notation adf.g = [/, 5 ] (Lie bracket of vector fields) and 
adhi.hj = {hi,hj} (Poisson bracket of Hamiltonian functions). 


Preliminary informal discussion. In order to compute singular controls, the usual method 
is to differentiate several times the switching function, until the control appears in a nontrivial 
way. If ||<l)(t)|| = 0 for every t G I, then hi{t) = h 2 {t) = 0, and, differentiating in t, we get, using 
the Poisson bracket, hi = {ho, hi} +M2{^2,^i} = 0 and /i2 = {ho, 112} + ui{hi,h2} = 0 along 
I. According to the Goh condition (see [16], see also below), if the singular arc is optimal, then 
the Goh condition {hi,h 2 } = {p, [gi,g 2 ]{x)) = 0 must be satisfied along I. Therefore we get that 
hi = {ho, hi} = (p, [f,gi]{x)) = 0 and h 2 = {ho,h 2 } = (p, [f,g 2 ]{x)) = 0 along I. 

Let us now assume that the vector fields gi and p 2 commute, i.e., [gi, 52 ] = 0. By differentiating 
again, we get 

hi = {ho, {ho, hi}} + Mi{hi, {ho, hi}} + U2{h2, {ho, hi}} = 0 , 
h 2 = {ho, {ho, h2}} + Mi{hi, {ho, h2}} + U2{h2, {ho, h2}} = 0 . 


If 


along I, then 


det Ai = det 


/{hi, {ho, hi}} 
\{hi, {ho, h2}} 


{h2, {ho, hi}}\ , 

{h2, {ho, h2}}y 


ui — ( — {ho, {ho, hi}}{h2, {ho, h2}} + {ho, {ho, h2}}{h2, {ho, hi}})/ det Ai, 
U2 = ({ho, {ho, hi}}{hi, {ho, h2}} — {ho, {ho, h2}}{hi, {ho, hi}})/ det Ai, 


( 15 ) 


and we say that the control u = {ui,U 2 ) is of order 1 (also called minimal order in [8, 12]). Note 
that ui and U 2 must moreover satisfy the constraint u\ + u\ ^ 1. Note also that, if moreover 
[ 51 , [ 7 , 52 ]] = 0 and [p 2 , [/,5i]] = 0, then (15) yields 


ui — —{ho, {ho, hi}}/{hi, {ho, hi}}, U2 — —{ho, {ho, h2}}/{h2, {ho, h2}}. 


Now, if {hi, {ho, hi}} = 0 and {h 2 , {ho, h 2 }} = 0 along I, then we must have {hi, {ho, hj}} = 0 , 
i,j = 1,2, i ^ j according to the Goh condition (see [16, 21], see also below), and hence we go on 
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differentiating. Assuming that [gi, [/, gi]] = 0 and [ 52 , [/iff 2 ]] = 0, we have 


[gi, adV-3i]] = bi, [/, ad/.gi]] 


-[/, [ad/.gi,gi]] - [ad/.^*, [gij]] = 0 , i = 1 , 2 , 


and we get 

= {^ 0 , ad^/iQ./ii} + U2{h2, ad^ho-hi} = 0, = {ho, ad^ho.h2} + ui{hi, ad^ho.h2} = 0. 

( 16 ) 

Due to higher-order necessary conditions for optimality (see below), an optimal singular control 
cannot appear in a nontrivial way with an odd number of derivatives, therefore we must have 
{/i 2 ,ad^/iQ./ii} = 0 and {hi,ad^/io./i 2 } = 0 along I. Accordingly, = 0 , i = 1 , 2 , gives the 
three additional constraints along the singular arc {ho, ad^ho-hi} = 0, {ho, ad^ho.h 2 } = 0, and 
{h 2 , ad^ho.hi} = —{hi, ad^ho.h 2 } = 0 . Derivating these constraints with respect to t, we get 

= ad'^ho-hi -|- iti{hi, ad^hg-hi} -|- tt 2 {ad^ho.hi, adho.h 2 } = 0 , 
h^"^^ = ad^ho.h2 -I- iti{ad^ho.h2, adho-hi} -|- U2{h2, ad^ho.h2} = 0 . 


Assuming that {h^, ad^ho-hi} < 0, i = 1,2 (generalized Legendre-Clebsch condition, see below) 
and that 


det A 2 = det 


/ {hi,ad^ho.hi} 
\{ad^ ho. h 2 , adho. hi} 


{ad^ho.hi, adho.h2} 
{h2, ad^ho.h2} 


^0 


along I, the singular control is given by 


ui = ( - (ad'^ho.hi){h 2 , ad^ho.h 2 } -b (ad^ho.h 2 ){h 2 , ad^ho-hi})/ det A 2 , 
U 2 = ((ad^ho.hi){hi, ad^ho.h 2 } — (ad^ho.h 2 ){hi, ad^hg-hi})/det A 2 . 


We say, then, that the singular control u = (mi,U 2 ) is of intrinsic order two. 


Precise definitions. Now, following [14], let us give a precise definition of the order of a singular 
control. 


Definition 1. 

of order q if 


The singular control u = {ui,U 2 ) defined on a subinterval I C [0,t/] is said to be 
d d'^ 


dui dt^ 
d 


(hi) =0, fc = 0,1, • • • ,2g - 1, 

d 


dut dt'^i 


{hi) ^ 0 , det 


du df 29 


4> ^0, f = l,2. 


along I. The control u is said to be of intrinsic order q if, moreover, the vector fields satisfy 


[gi,ad’"f.gi\ = 0, k = 1, ■ ■ ■ ,2q - 2, i = l,2. 


The condition of a nonzero determinant guarantees that the optimal control can be computed 
from the 2q-th time derivative of the switching function. Note that, in the definition, it is required 
that the two components of the control have the same order. 

We next recall the Goh and generalized Legendre-Clebsch conditions (see [16, 19, 21]). It is 
worth noting that in [ 21 ], the following higher-order necessary conditions are given even when the 
components of the control u have different orders. 


10 




Lemma 1. (higher-order necessary conditions) Assume that a singular control u = (mi, U 2 ) defined 
on I is of order q and is optimal. Then the Goh condition 

d df 

must be satisfied along I. Moreover, the matrix of which the {i,j)-th component is 

is symmetric and nonpositive along I (generalized Legendre-Clebsch Condition). 

In the problem (MTCP), as we will see, it happens that singular controls are of intrinsic 
order 2 , and that [ 171 , 52 ] = 0 , [ 51 , [ 7 , 52 ]] = 0 , and [ 52 J/, 5 i]] = 0 , so that the conditions given 
in the above definition yield [ 5 i,[/, 5 i]] = 0 , [ 52 , [ 7 , 52 ]] = 0 , [ 5 i,ad^ 7 . 5 i] = 0 , [ 52 ,ad^ 7 . 52 ] = 0 , 
(p, [ 51 , ad^ 7 . 51 ](a;)) 0 , (p, [ 52 , adV- 52 ](a;)) 0 , and 

(p, [ 5 i,ad^ 7 . 5 i](a;))(p, [ 52 , ad^ 7 . 52 ](a;)) - (p, [ 52 , ad^ 7 . 5 i](a;))(p, [ 51 , ad^ 7 . 52 ](a;)) y^ 0 , 

and we have the following higher-order necessary conditions, that will be used in the study of the 
problem (MTCP). 

Corollary 1. We assume that the optimal trajectory x{-) contains a singular arc, defined on the 
subinterval I of [0,tf], associated with a control u = {ui,U 2 ) of intrinsic order 2. If the vector 
fields satisfy [ 51 , 52 ] = 0, [gi, [f, gj]] = 0, for i,j = 1, 2, then the Goh condition 

(p(t), [ 5 i,ad 7 . 52 ](a;(t))) = 0 , (p(t), [gi,ad^f.g 2 ]{x(t))) = (p(t), [ 52 , adV- 5 i](a;(t))) = 0 , 

and the generalized Legendre-Clebsch condition (in short, GLCC ) 

{p{t),[g„ad^f.g^]{x{t))) 0, i = 1,2, 

(p(t), [ 5 i,ad^ 7 . 52 ](a;(t))) = (p(t), [ 52 , adV- 5 i](a;(t))) 

must be satisfied along I. Moreover, we say that the strengthened GLCC is satisfied if we have a 
strict inequality above, that is, {p(t), [gi,ad^f.gi\{x{t))) < 0. 

Corollary 1 follows from Lemma 1 and from the arguments developed in the previous informal 
discussion. It will be used in Section 4.2. 

We next investigate the singular junctions for the problem (13), and the chattering phenomenon. 

3.3 Chattering phenomenon 

One can find in [27] some results on the junction between an optimal regular arc and an 
optimal singular arc, for single-control affine systems, among which a result stating that, if the 
singular arc is of even order and if the control is discontinuous at the junction, then the junction 
must be nonanalytical (meaning that the control is not piecewise analytic in any neighborhood of 
the junction). In [30, 35], it is proved that such a nonanalytical junction between a regular arc 
and a singular arc of intrinsic order two causes chattering (see also [36]). When the control takes 
values in the unit disk, explicit analytic expressions for some optimal trajectories of linear-quadratic 
problems were given, e.g., in [24, 35]. However these results cannot be applied to (MTCP) because 
the control system is bi-input and the cost functional is the time; they are anyway a good source of 
inspiration to establish the results of that section. The following result is valid for general bi-input 
control-affine systems. 
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Theorem 1. Consider the optimal control problem (13). Let (x(-),p(-),p^,u(-)) be an optimal 
extremal lift on [0,1/]. We assume that u is singular of order two along an open interval I C [0,1/], 
and we denote this control by Us = (mis,U 2 s). We assume that j[ws(l)j| < 1 (i.e., the singular 
eontrol does not saturate the constraint) and that ■^^■^h 2 {x{t),p{t)) = 0 along I. Then the 
optimal control u must switch infinitely many times at the junction with the singular arc. In other 
words, there is a chattering phenomenon, which is due to the connection of a regular arc with a 
singular arc of higher-order. 

Proof. Since the singular control is of order two, it follows from the definition that 

d di^ d 

^^^i(a;(l),p(l)) = 0, fc = 0,---,3, * = 1,2, —-^h^{x{t),p{t))^0. 

Thus, we get from ■^^hi{x{f),p{f)) ^ 0 and Lemma 1 that 


3 

^^/ii(a;(l),p(l)) = 0, fc = 0. 


hj = 1,2, * j, 


and 


dm dC 


hr{x{t),p(t)) < 0, 




along the singular arc I. By assumption, we have ■^■^h 2 {x(t),p{t) = 0, and hence we can write 

hf\x{t),p{t)) = aioix{t),p{t)) Uisaii{x{t),p{t)) with aii{x{t),pit)) = ^^lij(x(l),p(l)) < 0. 

Without loss of generality, we consider a concatenation of a singular arc with a regular arc at 
time T G /. Assume that for some £ > 0 the control u is singular along (—£ + r, r), and that, along 
(r, T + £), the control u = (ui, 1 x 2 ) is given by m = /ii/jj<i>j] > 0, * = 1, 2. It can be easily seen from 
the assumption that jlrts]] < 1 that there exists at least one component of the singular control that 
is smaller than the same component of the regular control, i.e., Uks < Uk for fc = 1 or /c = 2. Then, 
it follows that 


= 0 ‘ko{x{T),p{T)) + Mfc(r)afcfc(x(r),p(r)) 

< ako(x(T),p(T)) + Uks(T~)akk(x(T),p(T)) = h^k\'^~) = 0 - 


(17) 


Hence the switching function hk has a local maximum at 1 = r and is nonpositive along the 
interval (t, t + £). It follows from the maximization property of the Hamiltonian that Uk ^ 0. 
This is a contradiction. If, instead, we assume m = /ii/ll<i>l] ^ 0, * = 1,2 over (r,r + £), then 
there must exists a control component Uks that is larger than Uks, i-e., Uks > Uk, and then we 
obtain h^^\T) > hl^\x~) = 0, which yields Uk ^ 0 and thus a contradiction. Then, if we assume 
Ui = < 0 and Uj = > 0, i,j = 1, 2, * /I j, we will have either Uig < m which gives 

a contradiction with the sign of m, or Uis ^ m and Ujg < Uj which gives a contradiction with the 
sign of Uj. A similar reasoning can be done for regular-singular type concatenations. 

Recall that the extremal is said singular if ll‘i?(t)j[ = \/hf{t) -I- h^t) = 0, t G I. Thus, the 
obtained contradiction indicates that the concatenation of a singular arc with a regular arc violates 
the PMP and thus there exists a chattering arc when trying to connect a regular arc with a singular 
arc. □ 


Remark 1. Note that, in this result, we have assumed that j]Msj| < 1. In the (nongeneric) case 
where the singular control saturates the constraint, in order to get the same result we need to 
assume that the strengthened GLCC is satisfied at the junction point, i.e., aii(a;(r),p(r)) < 0, and 
the control is discontinuous at the singular junction. 
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In addition, we have assumed that = 0. Actually, if -^■^h 2 {x{t),p{t)) ^ 

0, then singular and regular extremals can be connected without chattering. For example, (17) 
gives 

= «feo(a;(r),p(r)) + Ufc(r)afefc(a:(T),p(T)) + Mm(T)afcm(a;(T),p(r)) 

^ akoix{T),p{T)) + Uks{T~)akk{x{T),p{T)) + Um{T)akm{x{T),p{T)) 

= + akm{x{T),p(T)){Um{T) - Mms(r)) = Ofcm (^(t) , p(t)) («„ (t) - 

where akm(x(T),p(T)) = -^—■^hk{x{t),p{t))^ k,m = 1,2, k ^ m. In contrast to the previous 
reasoning, now the fact that akm(x{T),p(T)){um — Ums) > 0 does not raise any more a contradiction. 

In the next section, we analyze the regular, singular and chattering extremals for the problem 
(MTCP) by using the results presented previously. 

4 Geometric analysis of the extremals of (MTCP) 

In this section, we classify the switching points by their contact with the switching surface, and 
we establish that the optimal singular arcs of the (MTCP), if they exist, cause chattering. 

4.1 Regular extremals 

Normal extremals. Here, we consider normal extremals and we take p^ = —1. Let us consider 
the system (9), with the vector fields /, gi and g 2 defined by (10). Denoting the adjoint vector by 
P = {.Pv„,Pvy,Pv,,,Pe,Piij,P<i>,Puj„,Pujy), the adjoint equations given by the PMP are 

= 0 , Pvy = 0 , Pv^ = 0 , 

P 0 = -a cos tpipv,, cos 9 - py^ sin 9), 

p^ = a sin ij) sin 9py^ + a cos ipPvy + a cos 9 sin 'tjjPv,, — sin ip {uix sin (p + ojy cos (p) / cos^ ippg 

- {LOxS,i^(p + ^yCOS(p)/COS^ Ipp^, 

p^ = — {uJx cos (p — ujy sin (p) / cos ippg + {uJx sin (p + LUy cos (p)p.,f, 

— ta.nip[uix cos (p — uiy sin (p)p^, 

Pw* = — sinip/ cos ippe — cos (pp^j; — sin ip sincp/ cos ipp^, 

Pujy = — cos <p/ cos Ippg + sin pp^ — sin ip cos p/ cos ipp^j,, 

with the transversality condition 

Px^ (tf) sin 9 f cos Ip f - py^{tf)sinpf + py,^{tf) cos 9 f cosip f = 0. (19) 

The switching function is $(t) = {hi{t),h 2 {t)) = (bpyj^{t),—bp^^{t)) and is of class . The 
switching manifold P is the submanifold of of codimension two defined by 

r = {z = {x,p) e I Pyj,^ = Pojy = 0}. 

Let us fix an arbitrary reference regular extremal z(-) = {x{-),p{-)) of the problem (MTCP). 

If z{-) never meets P, then the extremal control is given by (14) along the whole extremal. 

If it meets P then there is a singularity to be analyzed. It is not even clear if the extremal 
flow is well defined when crossing such a point (we could lose uniqueness). Let us assume that the 
extremal z(-) meets P at some time to, and we set Zq = z{to) = (xo,po). Following [6, 7, 22], we 
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classify the regular extremals by their contact with the switching surface, i.e., if = 0 

and 7 ^ 0 for some k € N*, then the point zg is said to be a point of order k. Without loss 

of generality, we assume that to = 0- 

Let us analyze the singularity occuring at points of order 1, 2, 3 and 4 for the problem (MTCP). 


Points of order 1. 


Lemma 2. We assume that zg is of order 1. Then the reference extremal is well defined in a 
neighborhood of t = 0, in the sense that there exists a unique extremal associated with the control 
u = {ui,U 2 ) passing through the point zg. The control turns with an angle tt when passing through 
the switching surface L, and is locally given by 


Ul{t) 


ai t 

\/af + a| 1^1 


o(l), 


U 2 {t) 


02 t 
y^al + al 1^1 


0 ( 1 ), 


with ai = {hg,hi}{zg), 02 = {hg,h2}{zg). 

Proof. In the problem (MTCP), the vector fields gi and g 2 (defined by (10)) commute, i.e., 
[ 51 ,^ 2 ] = 0. By derivating the switching function, we get hi = {hg,hi} and ^2 = {hg,h 2 } along 
the extremal, and then, since the point is of order 1 , we have, locally, hi{t) = ait + oft) and 
/i 2 (t) = a 2 t + o(t), with af+ 02 ^ 0, and we also havepij^(t) = —a 2 t/b + o{t), Pujy{t) = ait/b + o{t), 
\/Pi^xW^ +Pi^y{t)^ = -\/af + al\t\/b + oft). The expression of the optimal control in the lemma 
follows. At the crossing point, both control components change their sign, i.e., Ui{0~) = —iti(0+), 
i = 1,2, which means that the control direction turns with an angle tt when crossing L. □ 


Points of order 2. 


Lemma 3. We assume that zg is of order 2. Then the reference extremal is well defined in a 
neighborhood oft = 0, in the sense that there exists a unique extremal associated with the control 
u = (ui,M 2 ) passing through the point zg. Moreover, the switching function is of class C°° in the 
neighborhood oft = 0, the control is of class C°°, and we have 


Ui{t) 


ai 

\/af+af 


0(1), 


U 2 ft) 


012 

\/a\ + al 


0(1), 


with ai = {hg,{hg,hi}}(zg) and 02 = {hg,{hg,h2}}{zg). 

Proof. The vector helds /, gi and g 2 , defined by (10), are such that [gi, [/, gj]] = 0, for i,j = 1,2 
(see Lemma 5 ). Then, according to the calculations done in Section 3.2, we have hi = {hg, {hg, hi}} 
and h 2 = {hg,{hg,h 2 }}, and hence the functions t 1 -^ hi{t), i = 1,2 are of class at 0. Locally, 
we have hi{t) = ^ait"^ + o{t^) and h 2 ff) = ^ 02 ^^ + o(t^). The expression of the optimal control 
follows and the control is continuous. 

Differentiating again the switching function, we have h^i\t) = ad^hg.hi + U 2 {h 2 ,a.d^hg.hi} 
and h^\t) = a.d^hg.h 2 + Ui{hi,ad^hg.h 2 }, because [ffi, [/,[/,ffi]]] = 0, for i = 1,2 (see Lemma 
5). Since the control is continuous, the switching function is at least of class at 0. Hence, 
locally we can write hi{t) = ^ait'^ + ^Pit^ + o(t^) and /i 2 (t) = ^a 2 t^ + + o(t^), where 

Pi = ad^hg.hi{zg) + otj! yj oifp a^fhj, ad^ hg.hi}{zg) for i,j = 1,2 and i j. We infer that the 
control is at least of class at 0, with ui(0) = l(a|/3i — 2 aia 2 P 2 )/{ai + and U 2 ( 0 ) = 

\{aiP 2 — ‘2aiOi2Pi)!{oil + ai)^/"^. We get the smoothness by an immediate induction argument: 
for k > 2, assuming that the switching function is of class and that the control is of class 
C^~^, then the {k + l)-th time derivative of the switching function can be written as = 
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u\ 'termi + term 2 , where ternii, i = 1,2 are terms involving time derivatives of the control of 
order lower than k — 2. Hence the switching function is of class since u is of class 

and hi = where all coefficients can be computed explicitly. The {k — l)-th 

time derivative of the control can be computed using the coefficients a^^p, p = 2, • • • , fc + 1, and 
hence the control is of class The result follows. □ 


Points of order 3. 


Lemma 4. We assume that zq is of order 3. Then the reference extremal is well defined in a 
neighborhood of t = 0 , in the sense that there exists a unique extremal associated with the control 
u = {ui,U 2 ) passing through the point zq. If bi = &d?hQ.hi{z[)) 0 and c = {/i 2 , ad^ho.hi}(zo) ^ 

0 , then the switching function is of class in the neighborhood of t = 0 and the control is 
discontinuous when passing through the switching surface T locally, and we have 


Ui{f) = 


A” 


+ o(l), t < 0, Ui{t) = 


PI 


T o(l), t > 0 , 


where, setting d = \J{—(? + b'f + 6 |) and Cij = —biC^ + bibj + b^, 

bj = i, 2 , 


-{-lybjcd + Cij + _ {-lybjcd + Cij 


bl + bl 


b\ + bl 


Proof. Using (10) and (16), we have = aWho-hi + U 2 {h 2 , ad^/iQ./ii} and h'^’ = ad^ho./i 2 — 
wi{ft, 2 ,ad^/iQ./ii} (see also Lemma 5). Assuming that we have locally hi{t) = \j3it^ + o(t^), we 


(3) 


infer that u^t) = 


0 ( 1 ). By substituting the control into the expression of we get 


/3i = 61 + c 




hi+hi 


and 132 = i >2 — c 


Pi 


The result follows by solving for t > 0 and t < 0. □ 


Remark 2. If c = {/i 2 , ad^ho.hi}(zo) = 0, then the switching function is of class at 0 and the 
control turns with an angle tt when passing through the switching surface T. 


Points of order 4. We assume that zq is of order 4. If {/ 12 ,ad^/io./ii}(zo) ^ 0, then 0 = 
hi = bl + U 2 C and 0 = h^ = 62 — uic, where bi = ad ho■hi{zo), c = {h 2 ,ad ho.hi}{zo). We 
have Ml = —bifc and M 2 = 62 /c. If moreover b^ + bl = c^, which indicates that the control 
Ui = aij \/ oi\p oi\ can be a regular control according to the value of hi{t) = |ait^ + o(t^), i = 1,2, 
at time 0, we get that M 2 /M 1 = 02/01 = — 61 / 62 , sign(ai) = 62 /c and sign(o 2 ) = — 61 /c. Then 

ui(t) = sign(c)^=| 2 = + 0(1), U2{t) = -sign(c)^=|i= + o(l). 

V + °2 V + °2 


If 


Oi = h[^\zo) = {6o,6i}(zo) + mi{6i,6i}(zo) + M2({62, 6 i}(zo) + {6i,c}(zo)) 

+MiM2{hl, c}(Zo) + M2{62, c}(Zo), 

02 = 6^'‘^(zo) = {6 o,62}(zo) + Mi({6i,62}(zo) - {6i,c}(zo)) + M2{62, 62}(zo) 

-Mi{6i,c}(zo) - MiM2{62,c}(Zo), 

then the extremal is well defined in a neighborhood of t = 0 and the control is continuous when 
passing the switching surface T. 
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If {h 2 ,a.<l^ho.hi}{zo) = 0, then {ho,&i}(zo) = {hoM}{zo) = 0, = 0, {h 2 , 6 i}(zo) = 

0 and {hi,bi}{zo) = {h 2 ,b 2 }{zo) (see the proof of Lemma 6 further), and we have h^\zo) = 
ui{hi,bi}{zo) and h^\zo) = U 2 {hi, bi}{zo)■ Assuming that we have locally ^(t) = + 

o(t"^) = i?oi^(cos(aln|t|), sin(aln|t|)) + o(t^) (we identify C = for convenience), with Rq > 0, we 
get u{t) = $(t)/||$(t)|| = and 

= i?o(4 + ia){i + ia){2 + ia){l + + o(l), 

which leads to Ro{a^ — 35a^ + 24) = {hi, fciK^o) and i?o(—lOa^ + 50a) = {hi,bi}{zo). It follows 
that Rq = (hi, 6 i}(zo)/(Q!^ — 35a^ + 24) with a G (1370/391, —871/614} if {hi,bi}{zo) < 0, and 
a G {578/1493,-5650/453} if {hi,6i}(zo) > 0. It is clear that the uniqueness of the extremal 
when crossing the point zq does not hold true anymore. The switching function $(t) converges to 
(0,0) when t —^ 0, while the control switches infinitely many times when t —>■ 0. Indeed, we will 
see further that this situation is related the chattering phenomenon. 

Abnormal extremals. Abnormal extremals correspond to = 0 in the PMP. We suspect 
the existence of optimal abnormal extremals in the problem (MTCP) for certain (nongeneric) 
terminal conditions. In the planar version of the problem (MTCP) studied in [36], if the optimal 
control switches at least two times then there is no abnormal minimizer. We expect that the 
same property is still true here. We are able to prove that the singular extremals of the problem 
(MTCP) are normal (see section 4.2), however, we are not able to establish a clear relationship 
between the number of switchings and the existence of abnormal minimizers as in [36]. Thus, in 
our numerical simulations further, we will assume that there is at least one normal extremal for 
problem (MTCP) and compute it. Note moreover that Lemmas 2, 3 and 4 are also valid for 
abnormal extremals. 


4.2 Singular and chattering extremals 

Let us compute the singular arcs of the problem (MTCP). According to [ 8 ], the singular 
trajectories are feedback invariants since they correspond to the singularities of the end-point 
mapping. This concept is related to the feedback group induced by the feedback transformation 
and the corresponding control systems are said to be feedback equivalent (see [ 8 , Section 4] for 
details). Recall that, equivalently, a trajectory x(-) associated with a control u is said to be 
singular if the differential of the end-point mapping is not of full rank. The end-point mapping 
E : R" X R X L°°(0, -boo; R) R" is defined as by £’(a:o, t/, u) = x{xo,tf ,u) where 1 1 —>• x{xq, t, u) 
is the trajectory solution of the control system associated to u such that x{xq^ 0 , u) = xq. 

Hence, we can replace the vector fields gi and g 2 with gi = and ^2 = Let us make 

precise the Lie bracket configuration of the control system associated with the vector fields /, 'gi 
and 'g 2 - 


Lemma 5. We have 


9i = 


duly ’ 


= [51,52] = 0, 


^ d d d 

ad/.oi = — cos 6/ cosib— + sin 07 -- — tan 0 cos 6-— = 0 , 
aO o^ip ocj) 

^ d d d 

a,df.g 2 = - sin cj)/ cos ' 0 — - cos 0 — - tan 0 sin 0 — = 0 , 

d& dip d(p 

d , „ , . d 


d 


aWf.gi = -^x-^ + a(cos 5 cos 0 -b sin5 sin 0 sin 0 )—-a(cos0sin5 — sin 0 cos5 sin 0 )——, 


' 50 


dvx 


dv. 
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„ d 

f.g 2 = + a{cos9sm(j) 

dtp 


— sin 0 cos ^ sin'0) —-a cos 0 cos 0—— 

dtij. dti,. 


— afsin 0 sin 0 + cos 0 cos 0 sin 0) ——, 

ov. 


[ 9 i,B,df.gi] = [gi,adf.g 2 ] = [g 2 ,adf.gi] = [g 2 ,adf.g 2 ] = 0, 


ad^f.gi 


d d d 

LO^ril/COSljj— - UJa:^2-^ + tan0Oi — 


d d d 

— aujy cos 0 sin 0 ^ -h aujy sin ip — - aujy cos ip cos 9 

OVx 


dVy 


dv, ’ 


ad^f.g2 


LOyfll/ COS0^ 


d d 

- uJx^ 2 TrT - uJx tan0rii — 


■ dip 


— auix cos iP! 


d 


d 


lgi,ad^f.gi] = Ig 2 ,ad^f.g 2 ] = 0, 


[gi,ad^f.g 2 ] 


-[ff2,ad^/.gi] 


A 


4 „ d 

ad^/.5i = —a(ujl + u}y){cos pcos9 + sin 0 sin0 sin 0)—-acos0 sin0(w^ + Wy)^— 

OVx OV',, 


d 

+ a{ujl +a;^)(cos0 * sin0 - cos0sin0sin0)--h {ojI + uJxOjI) 

" avz 


A 


ad'‘/.g2 


—a(w^ + a;^)(cos0sin0 — cos 0 sin 0 sin 0)—- h a cos p cos ip (ujI + ujy )—— 

" OVx OVy 


+ a{uJx + w^) (sin 0 sin 0 + cos 0 cos 0 sin ip) 


dvA^ 


[0i,ad^/.5i] 


,.„d -id d ■ I \ / I 9 

—acos0sin0—-h asin0—-acos0cos07-- [uJx sin0)/ cosip — 

OVx OVy OVz o9 

.8 8 

- UJx COS0— - UJx Sin0tan0 —, 
8ip Op 


[Si,ad^/.g2] 


d 

— 2aWy (cos 0 sin p — cos p sin ip sin 0) —- h 2aujy cos p cos 0 

OVx 


8vii 


8 


+ 2aujy (sin p sin 0 + cos p cos 0 sin 0) —-h (—w 

8vz 



8 


8 


[( 72 , ad'^/.gi] = — 2aWx(cos0cos0 + sin0sin0 sin 0)— - 2aujx cos 0 sin0—— 


8vx 


8Vy 


d d 

+ 2aujx (cos 0 sin 0 — cos 0 sin 0 sin 0) —-h (Sw^ + ) —, 

8vz 8p 
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[52,ad^/-52] 


= —acosijjsmO— -h asintb— -acosi/icos^—-(a;„ cos(j))/ cosi/j— 

OVx OVy OVz oO 

■ I , ^ 

+ (jjy Sin ffl-—— (jjy COS 0 tan w—, 

dip Off) 

where Jli = Ux cos (j) — ujy sin 0 and O 2 = uJx sin 0 + Wy cos 0. Moreover, we have 

dimSpan(gi,g2,ad/.5i,ad/.g2,ad^/.gi,ad^/.g2) = 6. 

Lemma 6. In the problem (MTCP), let us assume that {x[-),p{-),p^,u{-)) is a singular arc along 
the subinterval I, which is locally optimal in topology. Then we have u = (ui, U 2 ) = (0, 0) along 
I, and u is a singular control of intrinsic order two. Moreover, the extremal must be normal, i.e., 
p° ^ 0, and the GLCC 


a + gx sin 9 cos tp — gySintp + g^ cos 6 cos 0^0, (20) 

must hold along I. 

Proof. Using Lemma 5, we infer from $ = 0 and $ = 0 that 

{p, 9 i{x)) =Pu,y = 0, {p,g 2 {x)) = 0, 

{p,adf.gi{x)) = —pe COS0/COS0sin0 — tan0cos0 = 0, (21) 

{p, adf.g 2 ix)) = —pe sin 0/ cosip — cos cp — p^f, tan0 sin 0 = 0, 


and from $ = 0, that 

{p, ad^ f .gi{x)) = — iOxPrj, + a(cos0cos0 + sin0 sin 0sin0)p„^ 

— a(cos0sin0 — sin0cos0sin0)p„^ =0, 

{p,ad^ f.g 2 {x)) = + a(cos0sin0 — sin0cos0sin0)p^^ — a cos (p cos ippv^ 

— a(sin0sin0 + cos0cos0sin0)p„^ =0, 


( 22 ) 


along the interval I. Since dimSpan(gi, 32 , ad/.^i, ad/.(; 2 j ad^/.( 7 i, ad^/.g 2 ) = 6 , the six equations 
in (21)-(22) are independent constraints along the singular arc. Therefore, writing = 0, we 
get from Theorem 1 that 

{p, [h,a.d^ f-hix)]) =P 4 > = 0, (p, [p 2 ,ad^/.pi(x)]) = = 0, 

(p, ad^/.pi(a:)) =pgLOxIli/ cosip — p.,pu}x^2 + P<i>^x tan 0121 — Pv^auJy cos0 sin0 

+ px^ au}y sin ip — py^ awy cos ip cos 0 = 0, (23) 

(p, ad^ f.g 2 ix)) = -pgujyllif cos0 +p^Wa ;122 - P<i}0Jx tan0Ui - px^aoJx cos0 sin0 
+ px^ aujx sin ip — px^ auix cos ip cos 0 = 0, 


These four constraints are dependent: they reduce to two functionally independent constraints. 
Hence, with (21)-(22), we have 8 independent constraints along I. Now, we infer from (21)-(22)- 
(23) that P(j^ = Puiy = 0 and pe = Py. = P<^ = 0 along I. Derivating pg = 0 and p^ = 0, we 
get 

p„^ = tan0p„^, pxy = -tanip/ cos9px„, (24) 

and derivating again, that 0 = 0 = 0. It follows that Wx = ujy = 0. Using that H = 0 along any 
extremal, we get 

—p°COS0COS0 

0- + gx sin 0 cos 0 — py sin 0 + g^ cos 0 cos ip ’ 
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Substituting (24) and (25) into {p,ad^f.gi) and {p, f.g 2 ), we get 

(p, ad^/.pi) = (p, ad'‘/.p 2 ) = 0, (p, [pi, ad^/.p 2 ](a;)) = (p, [ 52 , ad^/.pi](a;)) = 0, 

(p, [gi,a(ff.gi]{x)) = (p, [p2, adV-52](a;)) =- 

COS cos 0 

To prove that u is of intrinsic order two, it suffices to prove that (p, [pi, ad^/.pi](a;)) ^ 0 along /. 
We prove it by contradiction. If (p, [p), ad^/.pi](a;)) = 0, then necessarily p„^ = 0 and this would 
lead to p„^ = p^^ = 0. It follows then from H = 0 that p° = 0. We have obtained that (p,p°) = 0, 
which is a contradiction. 

The fact that u = (iti, U 2 ) = (0, 0) simply follows from the fact that 

Ml = —{hi, ad^ho./iil/ad^/iQ.hi, U 2 = —{h2,ad^ho.h2}/ad'^ho.h2. 

Besides, if p° = 0, then p„^ = 0 and p^^ = Pv^ = 0, which leads to (p,p°) = 0 and thus to a 
contradiction as well. Therefore, p° < 0 (i.e., the singular arc is normal), and then (20) follows by 
applying the GLCC of Corollary I. □ 

We define the singular surface S, which is filled by singular extremals of the problem (MTCP), 
by 


5 * = 




UJoo =U}y = 0 , 


Pv. = 


Pe = Pijj = p<j> = Pui^ = Puiy = 0, Pv^ 
—p° cos 9 cos Ip 


: + Px sin 9 cos Ip — gy sin ip + gz cos 9 cos ip ’ 


Pvy 


tan 9py ^, 

— tani/)/cos0p„^ !•. (26) 


We will see, in the next section, that the solutions of the problem of order zero (defined in Section 
5.1.1) live in this singular surface S. 

The following result, establishing chattering for the problem (MTCP), is a consequence of 
Theorem 1, Lemma 5 and Lemma 6. 


Corollary 2. For the problem (MTCP), any optimal singular arc cannot be connected with a 
nontrivial bang arc. There is a chattering arc when trying to connect a regular arc with an optimal 
singular arc. More precisely, let u be an optimal control, solution o/(MTCP), and assume that u 
is singular on the sub-interval (^ 1 ,^ 2 ) C [0,t/] and is regular elsewhere. Ifti > 0 (resp., if t 2 <tf) 
then, for every e > 0, the control u switches an infinite number of times over the time interval 
[ti—e,ti] (resp., on [t 2 ,t 2 + e\). 

This result is important for solving the problem (MTCP) in practice. Indeed, when using 
numerical methods to solve the problem, the chattering control is an obstacle to convergence, 
especially when using an indirect approach (shooting). The existence of the chattering phenomenon 
in the problem (MTCP) explains well why the indirect methods may fail for certain terminal 
conditions. 

Note that, in the planar version of the problem (MTCP) studied in [36]), one can give suffi¬ 
cient conditions on the initial conditions under which the chattering phenomenon does not occur. 
Unfortunately, we are not able to derive such conditions in the general problem (MTCP). 


5 Numerical approaches 

In this section, we design two different numerical strategies for solving the problem (MTCP): 
one is based on combining indirect methods with numerical continuation, and the other is based on 
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a direct transcription approach. The first one may be successfully implemented when dealing with 
solutions without chattering arcs, and the second one is more appropriate to compute solutions 
involving chattering arcs. However, both approaches are difficult to initialize successfully because 
the problem (MTCP) is of quite high dimension, is highly nonlinear, and moreover, as a main 
reason, the system consists of fast (Euler angles and angular velocity) and slow (orbit velocity) 
dynamics at the same time. 

The occurrence of chattering arcs is an obstacle to convergence. Especially for indirect methods, 
the chattering phenomenon raises an important difficulty due to the numerical integration of the 
discontinuous Hamiltonian system. Direct transcription approaches provide a sub-optimal solution 
of the problem that has a finite number of switchings based on a (possibly rough) discretization. 
Actually, in case of chattering, we are also able to provide a sub-optimal solution with our indirect 
approach, by stopping the continuation before it would fail due to chattering. Though the sub- 
optimal solutions provided in this way may be “less optimal” compared with those given by a direct 
approach, in practice they can be computed in a much faster way and also much more accurately. 


5.1 Indirect method and numerical continuation 

The idea of this continuation procedure is to use the (easily computable) solution of a simpler 
problem, that we call herefter the problem of order zero, in order then to initialize an indirect 
method for the more complicated problem (MTCP). Then we are going to plug this simple, 
low-dimensional problem in higher dimension, and then come back to the initial problem by using 
appropriate continuations. 

This method actually gives an optimal solution with high accuracy. The problem of order zero 
defined below is used as the starting problem because the orbit movement is much slower compared 
with the attitude movement and it is easy to solve explicitly. As well, it is worth noting that the 
solution of the problem of order zero is contained in the singular surface S filled by the singular 
solutions for the problem (MTCP), defined by (26). 


5.1.1 Two auxiliary problems 

Problem of order zero. We dehne the problem of order zero, denoted by (OCPO), as a “sub¬ 
problem” of the complete problem (MTCP), in the sense that we consider only the orbit dynamics 
and that we assume that the attitude angles (Euler angles) can be driven to the target values in¬ 
stantaneously. Thus, the attitude angles are considered as control inputs in that simpler problem. 
Denoting the rocket axial symmetric axis as e and considering it as the control vector (which is 
consistent with the attitude angles 6, if), we formulate the problem as follows: 

V = ae + g, V{0) = Vo, V(.tf)//w, ||rZl|| = 1, mint/, 

where iZ; is a given vector that refers to the desired target velocity direction. This problem is easy 
to solve, and the solution is the following. 

Lemma 7. The optimal solution of (OCPO) is given by 


1 / kw — Vq 
a i tf 


tf = 


—02 + — doioa 


Pv = 


—P^ 




withk = {Vo,w) -b {g,w)tf, ai = a? - \\{g,w)w - gW^, 02 = 2{{Vo, w){g,w) - (Vq,^), and 03 = 

-||(Ho,h;)w - HolP- 
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Proof. The Hamiltonian is H = +py{ae + g), and we havep„ = 0, withp„ = and 

H = 0 along any extremal. It follows that py ^ 0 (indeed otherwise we would get also p^ = 0, and 
thus a contradiction). Hence there are no singular controls for this problem. The maximization 
condition of the PMP yields e* = p),/||p),||, and hence the optimal control is a constant vector. 
Moreover, according to the final condition Vitf)//w, the transversality condition is py T w, hence 
{e^,w) = 0, and using V{tf) = Hq + (ae + g)tf = kw we get that e** = _ gf it follows 

from the transversality condition that k = {Vo,w) + {pw)tf. The expression of t/ follows, using 
that ||e*||^ = 1. Using that 77 = 0, we get □ 

Since the vector e is expressed in the launch frame as {e)R = (sin0cos'0, — sin'0, cos 6 *sin'0)^, 
the Euler angles 0* G (— 7 r/ 2 , 7 r/ 2 ) and 0* G are given by 

0* = arctan(ei/e 3 ), 0* = —arcsin(e 2 ), (27) 

where e* is the i-th component of e*, for i = 1,2,3. 

Using the definition (26) of the singular surface 5, we check that the optimal solution of 
(OCPO) is contained in S with 6 = 9*, if = ip* and (p = (p* {cp* is any real number). Therefore, 
the relationship between (OCPO) and (MTCP) is the following. 

Lemma 8. The optimal solution of the problem (OCPO) actually corresponds to a singular solu¬ 
tion of (MTCP) with the terminal conditions given by 

Va:{0)=V^^, Vy{0)= Vy^ , (0) = , 

0(0) =r, 0(0) = 00 ,0(0) = 00 w,(0)=0, w,;(0)=0, 

Vz{tf)smipfVy{tf) cos9f cosipf = 0, Vz{tf)sin9f — Vx{tf)cos9f = 0, (29) 

9{tf)=0*, 0 (ty) = 00 , 0 (ty) = 00 uJz:{tf)=0, ujy{tf) = 0. (30) 

Due to this result, a natural idea of numerical continuation strategy consists of deforming 
continuously (step by step) the terminal conditions given in Lemma 8 , to the terminal conditions 
(11)-(12) of the problem (MTCP). 

However, because of the chattering phenomenon, we cannot make converge the shooting method 
in such a strategy. More precisely, when the terminal conditions are in the neighborhood of the 
singular surface S', the optimal extremals are likely to contain a singular arc (and thus chattering 
arcs). In that case, the shooting method will certainly fail due to the difficulty of numerical 
integration of discontinuous Hamiltonian system. Hence, we introduce hereafter an additional 
numerical trick and we define the following regularized problem, in which we modify the cost 
functional with a parameter 7 , so as to overcome the problem caused by chattering. 

Regularized problem. Let 7 > 0 be arbitrary. The regularized problem (OCPR).^ consists of 
minimizing the cost functional 


ftf 

Cj = tf-\-j {ul-\-ul)dt, (31) 

Jo 

for the bi-input control-affine system (9), under the control constraints —10 01, 7 = 1, 2, 

with terminal conditions (11)-(12). Note that, here, we replace the constraint -I- 0 1 (i.e., 

u takes its values in the unit Euclidean disk) with the constraint that u takes its values in the 
unit Euclidean square. The advantage, for this intermediate optimal control problem with the cost 
(31), is that the extremal controls are then continuous. 
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The Hamiltonian is 


Hj = {pj{x)) +ui{p,gi{x))+U 2 {p,g 2 {x)) + p° {1 + juf +-fuj), (32) 

and according to the PMP, the optimal controls are 

ui{t) = sat(-l,- 6 p<^^(t)/( 27 p°), 1), U 2 it) = sat(-l,&p<^Jt)/( 27 p°), 1), (33) 

where the saturation operator sat is defined by sat(— 1 ,/(t), 1 ) = —1 if f{t) ^ — 1 ; 1 if /(t) > 1 ; 
and f{t) if -1 ^ f{t) < 1 . 

As mentioned previously, one of the motivations for considering the intermediate problem 
(OCPR)^ is that the solution of (OCPO) is a singular trajectory of the full problem (MTCP), 
and hence, passing directly from (OCPO) to (MTCP) causes difficulties due to chattering (see 
Corollary 2). The following result shows that when we embed the solutions of (OCPO) into the 
problem (OCPR)^ , they are not singular. 

Lemma 9. An extremal of (OCPO) ean be embedded into the problem (OCPR)-,, by setting 

u(t) = (0,0), e{t) = 9*, ip{t)='ip*, (j){t) = (j)*, uJx{t) = 0, ujy{t) = 0, 

pe{t)=0, p^{t) = 0, p^{t) = 0, Pujx{t)=0, p,^y{t) = 0, 

where 9* and if* are given by (27), with terminal conditions given by (28) and (29)-(30). Moreover, 
it is not a singular extremal for the problem (OCPR).^ . The extremal eguations for (OCPR).y 
are the same than for (MTCP), as well as the transversality conditions. 

Proof. It is easy to verify that the embedded extremal is an extremal of the problem (OCPR).y 
and that the transversality conditions are the same. The control is computed from (33) which 
maximizes the Hamiltonian Hj, and we have Hj = 0 with = — 1. It follows from the PMP that 
the extremal equations are the same than for (MTCP). Then, for the problem (OCPR).^ , we 

have = 7 P°- Note that, in this case, the control Ui, i = 1,2 is singular if = 0. Hence 
there is no normal singular extremal for the problem (OCPR).^ . From Lemma 7, it is easy to see 
that ^ 0 and thus the extremals of (OCPO) are not singular extremals of (OCPR).y . □ 

5.1.2 Strategy for solving (MTCP) 

Continuation procedure. The ultimate objective is to compute the optimal solution of the 
problem (MTCP), starting from the explicit, simple to compute, solution of (OCPO). We proceed 
as follows: 

• First, according to Lemma 9, we embed the solution of (OCPO) into (OCPR).^ . For 
convenience, we still denote by (OCPO) the problem (OCPO) seen in high dimension. 

• Then, we pass from (OCPO) to (MTCP) by means of a numerical continuation procedure, 
involving three continuation parameters: the first two parameters Ai and A 2 are used to pass 
continuously from the optimal solution of (OCPO) to the optimal solution of the regularized 
problem (OCPR).^ , for some fixed 7 > 0 , and the third parameter A 3 is then used to pass 
to the optimal solution of (MTCP) (see Figure 3). 

The parameter Ai is used to act, by continuation, on the initial conditions, according to 

0(0) = 0*(1 — Ai) + 00'^ij V'(O) = ^”*(1 ~ ^ 1 ) + ^^(0) = ^^’*(1 ~ -^i) + 00^11 
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Figure 3: Continuation procedure. 


Wx( 0 ) — a;*(l — Ai) + Wy( 0 ) — oj*{1 — Ai) + w^pAi, 

where to* = uj* = 0, </>* = 0, and 9*, 'ip* are calculated through equation (27). 

Using the transversality condition (19) and the extremal equations = 0, Pv^ = 0 and 
'Pv^ = 0 , the unknown can be expressed in terms of and py^ as 

Pvy = {Pv^ sinOf cosipf +Pv^ cos6f cosipf)/shiipf, 

and hence the unknowns of the shooting problem are reduced to py^, py^, ? 5 e( 0 ), py.( 0 ), p<^( 0 ), 
Pi.cJx(0)i tf- The shooting function S\^ for the Ai-continuation is defined by 

Sx^ = {pujAtf), P<^y{tf): Pe{tf), Pij{tf), P4,{tf), 

Vz{tf) sinipf + Vy{tf) cos9f cosipf, Vz{tf) sin9f — Vx{tf) cos0/), 

where with p^ = —1 is calculated from (32) and ui and U 2 are given by (33). In fact, from 

Lemma 6 , we know that a singular extremal of problem (MTCP) must be normal, and since we 
are starting to solve the problem from a singular extremal, here we assume that p^ = —1. 

Note that we can use Sx-^ as shooting function thanks for (OCPR).^ . For problem (MTCP), 
if S'ai = 0 , then together with uJx{tf) = 0 and ujy{tf) = 0 , the final point {x{tf),p[tf)) of the 
extremal is then lying on the singular surface S defined by (26) and this will cause the fail of the 
shooting. However, for problem (OCPR).^ , even when x{tf) G S, the shooting problem can still 
be solved. 

Initializing with the solution of (OCPO), we can solve this shooting problem with Ai = 0, and 
we get a solution of (OCPR).^ with the terminal conditions (28)-(29) (the other states at t/ being 
free). Then, by continuation, we make Ai vary from 0 to 1, and in this way we get the solution 
of (OCPR).y for Ai = 1. With this solution, we can integrate extremal equations ( 8 ) and (18) 
to get the values of the state variable at tf. Then denote := 9{tf), ipe '■= 'fp{tf), (pe '■= (Pi^f), 

■= i^x{tf) and LOye ■■= iXyptf). 

In a second step, we use the continuation parameter A 2 to act on the final conditions, in order 
to make them pass from the values 0 e, ipe, (pe, ^xe and ujye, to the desired target values 0 /, ipf, 
<pf, oJxf and iXyf. The shooting function is 

Sx^ — (p^xipf^ (1 X2')UJxe X2U)x^i a^y(^/) (f A 2 )a^ye 

0{tf) - (1 - \2)9e - MOf: 'ip{tf) - (1 - A2)'0e “ A2'l/'/, (p{tf) “ (1 “ X2)(pe “ A2^/, 
Vz{tf) sinipf +Vy{tf) cos9f cosipf, v z{t f) sin9 f — Vx{t f) cos9f, Hj{tf)y 

Solving this problem by making vary X 2 from 0 to 1, we obtain the solution of (OCPR).., with the 
terminal conditions ( 11 )-( 12 ). 

Finally, in order to compute the solution of (MTCP), we use the continuation parameter A 3 
to pass from (OCPR).^ to (MTCP). We add the parameter A 3 to the Hamiltonian and to 
the cost functional (31) as follows: 

rtf 

C^ = tf+^ (“1 + “ 2)(1 “ ^ 3 ) dt, 

Jo 
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H{tf,X3) = (pj) + {p,9i)ui + {p,g 2 )u 2 +p° +p°j(uj + u'^)(l - A 3 ). 

Then, according to the PMP, the extremal controls are given by Ui = sat(—1,Miej 1), * = 1,2, 
where 

bPujy -bpuj^ 

^le — I / 1 ^2e — Z / 

-2^07(1 - A3) + h\ 3 ^pl^+pl^ -2^07(1 - A3) + h\ 3 ^pl^+pl^ 

The shooting function S\^ is defined as replacing with Hj{tf,X 3 ). The solution of 

(MTCP) is then obtained by making vary A 3 continuously from 0 to 1. 

Remark 3. Note that the above continuation procedure fails in case of chattering (see Corollary 
2), and thus cannot be successful for any possible choice of terminal conditions. In particular, if 
chattering occurs then the A 3 -continuation is expected to fail for some value A 3 = A| < 1. But in 
that case, with this value of A 3 , we have generated a sub-optimal solution of the problem (MTCP), 
which appears to be acceptable and very interesting for practice. Moreover, the overall procedure 
is very fast and accurate. Note that the resulting sub-optimal control is continuous. 

5.2 Direct method 

We now propose a direct approach for solving the problem (MTCP), where the control is 
approximated by a piecewise constant control over a given time subdivision. The solutions derived 
from such a method are therefore sub-optimal, in particular when the control is chattering (and 
in such a case the number of switchings is limited by the time step). Note that this approach is 
much more computationally demanding than the indirect one. 

Since the initialization of a direct method may also raise some difficulties, we propose the 
following strategy. The idea is to start from the solution of the problem (MTCP) with less 
terminal requirements, which is easy to obtain with a direct method, and then we introduce step 
by step the final conditions (12) of the problem (MTCP). We implement this direct approach 
with the software BOCOP and its batch optimization option (see [5]). 

• Step 1: we solve the problem (MTCP) with initial conditions (11) and final conditions 

ujy{tf) = 0, e{tf)=ef, V^{tf)sin0f-Vy:itf)cos9f = 0. 

These final conditions are the ones of the planar version of (MTCP) in which the motion 
of the spacecraft is 2D (see [36] for details). Numerical simulations show that, with such 
terminal conditions, the problem (MTCP) is easy and fast to solve by means of a direct 
method (a constant initial guess for the discretized variables suffices to ensure convergence). 

• Then, in Steps 2, 3, 4 and 5, we add successively (and step by step) the final conditions 

Vz(tf) sinipf + Vy{tf) cos9f cosipf = 0, ipitf) = ' 0 /, = 4'fi 0Jx(tf) = Wa,/, and for 

each new step we use the solution of the previous one as an initial guess. 

At the end of this process, we have obtained the solution of the full problem (MTCP). Note again 
that this direct approach is much slower than the indirect one, and that the resulting control has 
many numerical oscillations (see numerical results in Section 6.2). 
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6 Numerical results 


The structure of the rocket is presented in Figure 2 (b). We assume that the thrust / is 
flexible, i.e., it can turn ±6° in all directions, and its thrust is around Tatt = 1400 kN. The other 
thrusts are fixed with a total thrust Ttot = 1 X 10® kN. The rocket mass is 800 t, the length of the 
rocket is /i- = 50 m and its radius is = 2.5 m. Considering the rocket as a cylinder, we have 
Ix = ly = m{Zrl + /^)/12 and Iz = mrH^. The parameters a and b in (7) and (4) are therefore 
a = Ttotim Ri 12 and b = ~ 0.02. 

During the atmospheric ascent phase, the velocity of the rocket remains between several hun¬ 
dreds m/s and around 1000 m/s. Let v = +Vy+v1 be the modulus of the velocity, and let 
ipy and 4>v be the flight path angles that we use to calculate the components of the velocity in Sn 
frame, i.e., = v sin Oy cos ipy, Vy = —vsimpy and Vz = v cos Oy cos ij^y. In this section, the initial 

values of the angles 9y and ipy are chosen equal to the initial values of the angles 9 and ip. This 
means that, before the maneuver, the rocket is on a trajectory with angle of attack equal to zero. 

In the numerical simulations, we set Vxg = vq sin 9o cos ipo , Vy^ = —vq sin ipQ , Vzg = vq cos 9^ cos ipo 
and take the other values needed in the initial condition (11) and the final condition (12) in the 
following table 


(TCI): 

e 

0 

II 

e 

0 

II 

0 

00 

= 

75°, 

Ipo = 

0.5°, 

(po = 

= 0° 


f ~ ^Vf ~ 

= 

85' 

", 

= 5° 

, <^/ = 

= 0° 


(TC2): 

e 

0 

II 

e 

0 

II 

0 

00 

= 

70°, 

'ipo = 

0.5°, 

(po = 

= 0° 


f ~ ^V.f ~ 

= 

85' 

", V'.f 

= 5° 

, 

= 0° 


(TC3): 

o' 

II 

0 

3 

II 

0 

3 

00 

= 

85°, 

■00 = 

0.5°, 

(po = 

= 0° 

^X 

0 " 

II 

3 

II 

= 

75' 

", 

= 5° 

, 

= 0° 



Table 1: Terminal conditions 


Note that ?;o is the module of velocity at time 0. In the next two subsections, we will choose 
different values of Vg, and so here we do not assign to it a specific value. Moreover, we set 7 = 50 
as the weight of the L 2 -norm control term in the cost functional of problem (OCPR).^ . 

6.1 Numerical results without chattering 

The indirect method combined with numerical continuation described in Section 5.1 is imple¬ 
mented using a predictor-corrector continuation method, where the prediction is made thanks to a 
Lagrange polynomial. The Fortran routines hybrd.f (see [28]) and dop853.f (see [18]) are used, 
respectively, for solving the shooting problem (Newton method) and for integrating the ordinary 
differential equations (with prediction). 

The Euler angle 9 is usually called the pitch angle., and a pitching up maneuver designates a 
maneuver with terminal condition 9f > 9g, while a pitching down maneuver designates a maneuver 
with terminal condition 9f < 9g. 

Pitching up maneuvers. We set vg = 1000 m/s and we use the numerical values denoted 
by (TCI) in Table 1. The components of the state variable are reported on Figure 4. The 
optimal control, the adjoint variables it) and p^y^ (t) and the modulus of the switching function 
<I>(t) = b{piy ^, —pu !^) are reported on Figure 5. We observe that the optimal control switches twice, 
at times 8.8 s and 25.8 s. These two switching points are of order 1 (i.e., <I>(t) = 0 and $(t) 0). 

Accordingly with Lemma 2, the control turns with an angle n at those points. 
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Figure 4: State variable and the optimal control with (TCI) and = 1000. 




Figure 5: Adjoint variable and the switching function with (TCI) and vq = 1000. 


Let us give another numerical example, taking the same terminal conditions as previously 
except for uq, and we take Vg = 1500 m/s. The time history of the state, of the optimal control 
and of the switching function are reported on Figures 6 and 7. One can see on Figure 7 that the 
optimal control turns two more times with an angle tt due to two switching points of order one. 

Pitching down maneuvers. We set vg = 1500 m/s and we use the numerical values denoted 
by TC3 in Table 1. The optimal solution is drawn on Figures 8 and 9. The shorter maneuver 
time tf indicates that it is easier to turn clockwise the axis of the velocity vector than to turn it 
anti-clockwise. This corresponds to the intuition. The reason is that the total force induced by 
the gravity force tends to reduce Vx, i.e., it helps the velocity to turn clockwise, and so together 
with the rocket thrust force, the maneuver time is less than that of the anti-clockwise case. 

Note that the derived time history of the adjoint variable (for both pitching up and pitching 
down maneuvers) do not have the same order of magnitude, i.e., Poj^ and are ten times larger 
than pff and p^, and are thousand times larger than and py^. This indicates again that 

the shooting method is difficult to initialize successfully. 

We note that the indirect strategy proposed in Section 5.1.2 is efficient also because the smallest 
adjoint variables Pv ^, py^ and py^ are already quite accurately estimated thanks to the problem of 
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Figure 6: Time histories of state with (TCI) and uq = 1500 m/s. 




Figure 7: Time histories of control and switching function with (TCI) and vq = 1500m/s. 


order zero. 

6.2 Numerical results with chattering arcs 

Sub-optimal solution by the indirect approach. On Figures 10 and 11 is given a sub-optimal 
solution of (MTCP) with the terminal conditions (TC2) of Table 1 and vq = 2000 m/s. Due to 
chattering, the continuation parameter A 3 stops at value Ag = 0.98 (see Remark 3). Observing 
from Figure 11, the switching function pass four times the switching surface F are small between 
time 26.5 and 40.3. The control, instead of bang-bang or singular, is continuous. The cost of this 
trajectory is 69.3 and the final time tf = 66.0 s. 

Sub-optimal solution by the direct approach. With the same terminal conditions as above, 
we now use the direct method described in Section 5.2. Numerical simulations show that the 
initialization step for the direct method procedure is quite robust (a constant initial guess is 
enough). The results are reported on Figures 12, 13 and 14. 

We observe that, when t G [23,43], the control oscillates much with a modulus less than 1: 
this indicates that there is a singular arc in the “true” optimal trajectory, and therefore chattering 
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Figure 8: Time histories of state with TC3 and vq = 1500. 




Figure 9: Time histories of control and switching function with TC3 and vq = 1500. 


according to Corollary 2. 

Note that, along the singular arc, the variables uJx, ujy, , pg, p^ and p^ are almost equal 

to 0, and we check that this arc indeed lives on the singular surface S defined by (26). Therefore, it 
turns out that there is a singular arc in the optimal trajectory, causing chattering at the junction 
with regular arcs. 

The maneuver time is tf = 65.4 s. Compared with that of the sub-optimal solution derived 
from the indirect strategy, only 0.6 s are gained with the direct method. The direct approach 
is hundreds of times slower than the indirect approach and the obtained control presents many 
oscillations, which is not much appropriate for a practical use. 

On Figure 4, 6, 8, 10 and 12, we note that the attitude angles first tend to reach the values 
0* deg, '0* (i.e., 9* = 176.9 deg and ip* = 18.5 deg for Figures 4 and 6; 9* = —17.4 deg and 
Ip* = 25.2 deg for Figure 8; 0* = 176.1 deg and ip* = 11.2 deg for Figures 10 and 12), and then 
turn back to reach their final values. Actually, doing more numerical simulations with different 
terminal conditions (note reported here), we observe that the extremals have a trend to first go 
towards the singular surface and then to get back to the target submanifold. We suspect that this 
is due to a turnpike phenomenon as described in [34], at least when the required transfer time is 
quite large. 


28 



















































Figure 10: Time histories of state with (TC2) and vq = 2000. 




Figure 11: Time histories of control and switching function with (TC2) and vq = 2000. 


7 Conclusion 

We have studied the time optimal control of the rocket attitude motion combined with the 
orbit dynamics. The problem (MTCP) is of interest because of the coupling of guidance and 
navigation systems. However, this problem is difficult to solve because of the occurence of the 
chattering phenomenon for certain terminal conditions. 

Using geometric control, we have established a chattering result for bi-input control-affine 
systems. We have also classified the switching points for the extremals of the problem (MTCP), 
according to the order of vanishing of the switching function, showing the behavior of the control 
at the singularities. 

In order to compute numerically the solutions of problem (MTCP), we have implemented 
two approaches. The indirect approach, combining shooting and numerical continuation, is time- 
efficient when the solution does not contain any singular arcs. For certain terminal conditions, 
the optimal solution of (MTCP) involves a singular arc that is of order two, and the connection 
with regular arcs can only be done by means of chattering. The occurrence of chattering causes 
the failure of the indirect approach. For such cases, we have proposed two possible numerical 
alternatives. Since our indirect approach involves three continuations, one of them being concerned 
with a continuation on the cost function (and thus on the Hamiltonian and the control), we have 
proposed, as a first alternative, to stop this last continuation before its failure: in such a way, we 
obtain a sub-optimal solution, which seems to be very acceptable for a practical use. The second 
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Figure 12: State variable x{t) (vg = 2000 m/s). 



Figure 13: Optimal control and ||'b(t)|| (uq = 2000m/s). 


alternative is based on a direct approach, and then we obtain as well a sub-optimal solution having 
a finite number of switchings, this finite number being limited by the chosen step of the subdivision 
in the discretization scheme. In any case, the direct strategy is much more time consuming than 
the indirect approach. Note that, in both cases, it is not required to know a priori the structure 
of the optimal solution (in particular, the number of switchings). 

As an open issue, one may consider to add atmospheric forces in the model. Since the magnitude 
of the aero-forces is low (at least, it should be much smaller than the rocket thrust), we expect 
this extension to be doable, for instance by means of an additional continuation. 
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